Wind and boundary layers in Rayleigh-Benard convection. 
I. Analysis and modeling 
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The aim of this paper is to contribute to the understanding and to model the processes controlling 
the amplitude of the wind of Rayleigh-Benard convection. We analyze results from direct simulation 
of an L/H = 4 aspect-ratio domain with periodic sidewalls at Ra = {10 s , 10 6 , 10 7 , 10 s } and at 
Pr = 1 by decomposing independent realizations into wind and fluctuations. It is shown that deep 
inside the thermal boundary layer, horizontal heat-fluxes exceed the average vertical heat-flux by 
a factor 3 due to the interaction between the wind and the mean temperature field. These large 
horizontal heat-fluxes are responsible for spatial temperature differences that drive the wind by 
creating pressure gradients. The wall fluxes and turbulent mixing in the bulk provide damping. 
Using the DNS results to parameterise the unclosed terms, a simple model capturing the essential 
processes governing the wind structure is derived. The model consists of two coupled differential 
equations for wind velocity and temperature amplitude. The equations indicate that the formation of 
a wind structure is inevitable due to the positive feedback resulting from the interaction between the 
wind and temperature field. Furthermore, the wind velocity is largely determined by the turbulence 
in the bulk rather than by the wall-shear stress. The model reproduces the Ra dependence of wind 
Reynolds number and temperature amplitude. 

PACS numbers: 44.25.+I, 47.27.ek, 47.27eb, 47.27.te 
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I. INTRODUCTION 



One of the characteristic features of Rayleigh-Benard 
convection is a large scale circulation or 'wind, which is 
generated autonomously by the system and is of great 
importance for the effectivity of the heat transfer [lj . Al- 
though first observed in a large aspect-ratio T = L/H cell 
0, the wind has been studied mostly in smaller aspect- 
ratio cells @, i, i, H 0, H ©, EI HO, El ■ The wind has 
complex dynamics, in that it changes its direction errat- 
ically at timescales far exceeding the convective turnover 
time [1, 0. In the case of cylindrical cells, there are 
two separate ways for reversals to occur [l3|, EH- First, 
the wind structure can change its orientation by rotat- 
ing in the azimuthal direction, which leads to reversals 
if the system rotates over 180°. The second mechanism 
for reorientation is by cessation, when the large scale cir- 
culation briefly halts and restarts with a different ran- 
dom orientation. The wind dynamics change depending 
on the aspect-ratio T and the Raylcigh number Ra. In 
cylindrical T — 1/2 domains, the wind structure (nor- 
mally one roll throughout the entire domain) breaks up 
into two counter-rotating rolls on top of each other [15[ 
around Ra = 10 10 . At even higher Ra, roughly around 
10 12 , the wind substantially weakens (2, [lE El- For 
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large aspect-ratio domains, the wind structure tends to 
be weaker relative to the fluctuations [H, El, HI HH HI • 

Several models have been developed recently to ex- 
plain the complex long-term dynamics of the wind, in 
particular the wind reversals and reorientations. The 
first model to explain wind-reversals was by Sreenivasan 
et al. 0, which is based on the conceptual image of a 
double-well potential representing the preference for an 
average clock- wise or counter-clockwise motion. The tur- 
bulence is modelled by stochastic fluctuations, which are 
responsible for sudden reversals when strong enough to 
overcome the energy barrier separating the two states. A 
different approach was taken by Fontenele Araujo et al. 
[23l | , who derived a deterministic model describing the dy- 
namics of a thermal on a circular trajectory in a linearly 
unstably stratified fluid. The resulting equations are sim- 
ilar to the Lorenz equations and exhibit chaotic flow re- 
versals in a specific region of the Ra-Pr phase space. The 
two-dimensional models described above can only repro- 
duce reversals by cessations, and do not facilitate reori- 
entation by rotations, which occur more often [l3|, EH- 
Brown and Ahlers [2J] recently presented a model which 
is capable of predicting reorientations both by rotations 
and cessations. This model is inspired by the Navier- 
Stokes equations and constitutes two stochastical differ- 
ential equations, one for the temperature amplitude and 
one for the azimuthal orientation. 

Despite these significant advances in the understand- 
ing of the long-term wind-dynamics, it is currently not 
clear exactly how the wind is driven and how the tur- 
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bulence and wall-fluxes influence the wind amplitude. It 
is known that the wind is sustained by the spatial dif- 
ferences in mean temperature along the sidcwalls (25[. 
However, it is not clear what generates these tempera- 
ture differences, and what the relation between the tem- 
perature differences and the wind velocity is. In this 
paper, we use direct numerical simulation of a rectangu- 
lar T = 4 domain at Pr = 1 and Ra = {10 5 , 10 6 , 10 7 , 10 8 } 
with periodic lateral boundary conditions to provide in- 
sight into these questions. We derive a model for the wind 
based on the Reynolds- averaged Navier-Stokes equations, 
which consists of two coupled ordinary differential equa- 
tions for the average wind velocity and temperature am- 
plitude. This simple conceptual model provides insight 
in the role of turbulence in the bulk on the wind velocity 
and the neccessity for a wind-structure to develop. In the 
accompanying paper (26j , we will focus on the boundary 
layers at the top and bottom walls and their interaction 
with the wind, and propose new scaling relations for X u 
and Cf. 

The paper is organized as follows. The governing 
equations, averaging strategies and their relation to 
the system's symmetries are discussed in section III Al 
The method of wind extraction by symmetry-accounting 
ensemble-averaging is outlined in section III Bt Similar 
to domains with sidewalls, a wind structure develops for 
unconfined domains [H, [2l|, [27|, H, Hg|. As the wind 
structure is not kept in place by side walls, it can be lo- 
cated anywhere in the domain because which complicates 
extracting the wind structure. However, by identifying 
the wind structure and proper alignment of different re- 
alizations (by accounting for symmetries), a wind struc- 
ture can also be unambiguously defined for unbounded 
domains [29[ . Details about the code and simulations are 
discussed in section III CI Some results of Nu and Re as 
a function of Ra are presented in section IIIII The wind 
and the temperature field following from the symmetry- 
accounted averaging are presented in section IIVI The 
decomposed profiles of kinetic energy are presented in 
section HVBI eliciting the importance of the wind for the 
dynamics of the flow. It turns out that the wind struc- 
ture has a significant influence on the redistribution of 
heat in the system, as is discussed in section IIV CI Sec- 
tion IIVDI contains a discussion how the wind is main- 
tained by a study of the momentum and temperature 
budgets at several positions of the flow, and a detailed 
feedback mechanism is sketched. Then, the findings are 
synthesized in a simple conceptual model in section IVl 
and conclusions are presented in section [VTl 



II. BACKGROUND 
A. Theory 

Raylcigh-Bcnard convection is generated when a layer 
of fluid with thickness H between two parallel plates is 
subjected to a positive temperature difference AO be- 



tween top and bottom plate. The positive tempera- 
ture difference causes the buoyant fluid to become unsta- 
ble, causing convection and thereby enhancing the heat- 
transport through the layer. In the dynamics one can 
observe organized motion such as plumes, jets and wind 
fl2| . For an incompressible Boussinesq fluid with isobaric 
thermal expansion coefficient (3, viscosity v and thermal 
diffusivity k, the governing equations are 

dtu, + OjUjU, = -p^ 1 d l p + vdjiit + (3g<d5 l3 , (1) 
d t Q + djUjQ = k<9?6, (2) 
djUj = 0. (3) 

Here p is the density, g the gravitational constant, itj rep- 
resents the fluid velocity, O the temperature and p the 
pressure. No-slip velocity and fixed temperature are en- 
forced on the top- and bottom walls. The problem can be 
characterized by the Prandtl number Pr = vk~ x which 
represents the ratio of viscosity and thermal diffusivity 
and the Rayleigh number Ra = /3gA9i/ 3 (^K) _1 which 
relates the buoyant and viscous forces. The system re- 
acts by convective motion characterized by the Reynolds 
number Re = UHv~ l and by an enhanced heat transfer 
through the Nusselt number Nu = (I)H(kAQ)^ 1 which is 
the non-dimensional heat-flux through the top and bot- 
tom wall. Here U is a characteristic velocity and <fi the 
heat-flux. Both Re and Nu are unknown a priori. 

Since definitions for the processes occurring in 
Rayleigh-Bcnard convection arc not unambiguous, a 
small glossary is given here. We prefer to use the term 
wind structure, which generalizes the terms wind and 
large scale circulation, in that it involves both the ve- 
locity and the temperature field. This wind structure 
normally features convection rolls, which are the quasi- 
steady roll-like structures. Thermals and plumes are the 
unsteady structures erupting from the boundary layers 
and propagating into the bulk. Spatial averages will be 
denoted by ()a and ()h for volume-, plane- and 

height-averaging, respectively. The plane-average is in 
the homogeneous (x and y) directions. Time and ensem- 
ble averages will be denoted by ()t and (). 

In what follows a domain of size LxLxH with L = TH 
and r the aspect-ratio will be considered. Periodic 
boundary conditions are imposed on the side walls. Ap- 
plying (),4 to the incompressibility constraint §5§ and us- 
ing impermeability at the top and bottom wall yields that 
the plane-averaged velocities (u)a = (v)a = (w)a = 0. 
Taking the ensemble average of the temperature equation 
([2]) and the fixed temperature boundary conditions gives 
after some manipulation that 

Nu=^(( w 'e')-«5,(e)) > (4) 

which states that the mean total heat-flux is constant in 
the vertical and directly related to Nu. 

Interesting differences exist in the standard way of 
averaging between experiments, simulation and theory. 
We focus on laterally unbounded domains or domains 
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with periodic boundary condition and will use the over- 
bar X to denote a generic averaging operator. Experi- 
ments normally employ the time-average (X) t and the- 
ory the ensemble average (X). In simulations of un- 
bounded Rayleigh-Bcnard convection it is customary to 
use a plane-average (X)a, because it can be evaluated 
at every time instant. The underlying assumption is 
that X coincides with the ensemble average (X) and the 
time average (X) t , but there are some subtleties that 
require attention here. It can be imagined that {X)a 
will approach (X) for T sufficiently large, as a typical 
realization is expected to be of size O(H) by which the 
domain would contain roughly T 2 of those realizations. 
The time average (X) t produces one independent real- 
ization every 0(t*) with t* = H/U the typical timescale, 
and it can be expected that for averaging over sufficiently 
long times it converges to the ensemble average so that 
(X) t = (X)a = (X). However, this presumes that the 
system's phase space is not partitioned, i.e that the sys- 
tem will visit all its possible states within finite time. 
When this condition is satisfied the system is ergodic, 
and this is one of the primary assumptions underlying 
turbulence theory [30L l3l| From the continuity equation, 
it follows that {ui) a = 0, by which all natural averages, 
i.e. long-time, ensemble and spatial averages vanish as 
u = v = w = 0. Hence one would conclude that Rayleigh- 
Benard convection is comprised purely of fluctuations, 
which is in conflict with the ubiquitous large scale circu- 
lation or wind. 

The paradox of the existence of a mean wind and the 
restriction o{u = v = w = can be resolved by taking 
into account the symmetries of the problem [23]. When 
there are symmetries in the domain, there is a chance 
for symmetric conjugate modes (such as clockwise and 
counter-clockwise mean flow in the cell) to cancel each 
other, given enough time (through wind-reversals) or re- 
alizations. By accounting for symmetries before perform- 
ing ensemble-averaging, all fields are properly 'aligned' 
before the averaging takes place, allowing the modes that 
would normally be cancelled by their symmetric conju- 
gates to persist. The resulting average field of velocity 
and temperature is the wind structure and in the fluctu- 
ations are the actions of the plumes. 



B. Symmetry-accounted ensemble-averaging 

The rationale of symmetry-accounted ensemble- 
averaging has been presented for general domains else- 
where [29[ and we discuss here only the application to our 
case with periodic side walls. The system has two sym- 
metries: a discrete rotational symmetry and a continuous 
translational invariance. The most important symmetry 
to take into account here is the translational invariance 
in x,y. When considering an ensemble of realizations 
{XW,XW, . . .,1^}, it can be expected that a wind 
structure is present in all of them, although its location 
will differ per realisation. When one takes the average 



of this ensemble, the wind structure will be averaged out 
so that nothing but fluctuations remain (Fig.[Th). How- 
ever, due to the translational invariance, one can trans- 
late a realization and obtain another valid solution to 
the equations. By translating each realization X^ over 
a distance d^ a > such that the wind structures become 
aligned, the averaging out of the wind can be prevented, 
as is sketched in Fig. [T|d. 

The translational operator can be denoted by Sd with 
d = (d x ,dy) representing the relative displacement. Op- 
erating on a field X, the translational operation is simply 
SdX = X(x — d x ,y — d y , z). Symmetry-accounted aver- 
aging then, means to translate each realization a before 
averaging as 

N 

" =1 (5) 

N y ' 

= Y / X ia) (x-d x a} ,y-d y a \z), 

a=l 

where is chosen such that the wind structure does 
not average out. An alternative way to look at symmetry- 
accounted ensemble-averaging is that it involves a prepro- 
cessing step before performing the ensemble-averaging. 
The fluctuating field is defined as 

X'W = X^ {x - d&*> , y - 4 Q > ,z)-X{x,y,z), (6) 

and it is straightforward to prove that X' = 0. Hence, 
the results can be interpreted exactly the same way as 
those from classical Reynolds-decomposition. 

The symmetry-accounted ensemble average X is 
closely related to the classical (ensemble, long-time or 
spatial) average X, and we will point out some useful 
relations between the two. Due to translation invari- 
ance all statistics X are a function of z only, whereas 
the symmetry-accounted average X retains the full three- 
dimensional structure. The first important relation is 
that the plane-average of the symmetry-accounted aver- 
age is identical to the classical average as 

(X)a = X (7) 

which follows directly from substitution of the two dif- 
ferent decompositions X — X(x,y,z) + X'(x, y, z) and 
X = X(z) + X"(x,y,z) into the expression (X)a- The 
second useful relation pertains the variance, and is given 
by 

{XX) a + (x^X 1 ) a = X~X~ + X ir X" (8) 

which can be obtained similarly. Expression |8]) is par- 
ticularly useful for the analysis of the profiles of kinetic 
energy pV B[) and for the decomposed vertical heat-fluxes 
(section HVCl) . 

If the wind structure was known a priori, the displace- 
ment d would be the only unknown per realization, and 
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FIG. 1: Ensemble averaging in domains with periodic side walls, a) Classical averaging results in zero mean wind; b) When 
accounting for symmetries by translating the realizations if necessary, the wind structure is preserved. 



([5]) could be applied immediately. Unfortunately this 
is not the case, as both the wind structure and d are 
unknown. Therefore, an iterative technique is used by 
which the wind structure and the displacements are de- 
termined simultaneously, gradually improving the esti- 
mation for the wind structure in successive iterations 
f32| . The only assumption needed for this method is 
that - among the majority of the realizations - only one 
persistent structure (mode) is present inside the domain. 

To start the iterative process a reference field -Xo(x) is 
needed, for which an arbitrarily picked realization is used 
- the wind structure is present in every realization so 
the starting point should not make a difference. Using a 
cross-correlation function C(X,Y), every realization can 
be compared to Xq(x), and the location of maximum 
correlation is picked as the displacement vector: 



tcrmined by 



d (Q) <- max C{S r X^ a \X ) 



(9) 



There is some freedom in choosing how to calculate the 
overall 2D (in x and y) correlation field, as it can be con- 
structed from any combination of the three-dimensional 
fields X <E {ui,<d,p}. In this case we opted for the in- 
stantaneous height-averaged temperature (&)h which is 
closely related to the wind structure as (&)h > where 
w > and vice versa. Denoting the reference field by 
Xq(x,u) = (&q)h and a different realisation by Y{x,y), 
the cross-correlation function is given by 



C(S r Y,X ) 



II Y '( x ~ r x,y- fy)X' Q (x, y)dxdy 



<tx&y 



(10) 



Here, X' = X - (X Q ) A and Y' = Y-(Y) A are the devi- 
ations from the mean, and ax and ay are the standard 
deviations of Xq and Y. The displacement vector d is 
just the coordinate pair (r Xl r y ) for which the correlation 
is maximal. For computational efficiency, the correlation 
is determined via FFT's. After calculating d( Q > for all 
realizations, a new and improved estimation can be de- 
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(ii) 



Repeatedly applying ([9]) and (fTTj) with Xq replaced with 
X n and until X n+ \ = X„ = X results in the wind struc- 
ture, or symmetry accounted average, as well as the rela- 
tive displacements d^"). It is emphasized that vertically 
averaged fields are only used to determine the relative 
displacements d^"); the resulting wind structure is fully 
three-dimensional. 



C. Simulation details 

Direct numerical simulation (DNS) is used to generate 
the independent realizations for the symmetry-accounted 
averaging. The code is based on finite volumes and has 
the equations ([T][3]) discretized and implemented on a 
staggered grid. Central differences are used for the spa- 
tial derivatives and time integration is by a second order 
Adams-Bashforth scheme. The code is fully parallelized 
and supports grid clustering in the wall-normal direc- 
tion. Special attention has been given to conservation 
of variance by preserving the symmetry- properties of the 
discrete advective and diffusive operators 
details of the code can be found elsewhere 

Resolving all the length-scales makes direct numeri- 
cal simulation a powerful research tool, as one has the 
complete four-dimensional solution of the Navier-Stokcs 
equations at hand. However, DNS is limited to relatively 
low Re as the computational demands quickly become 
prohibitive, scaling approximately as Re 3 . Furthermore, 
both the thermal and hydrodynamic boundary layer, Ae 
and X u respectively, should be fully resolved as under- 
sampling will lead to overestimation of Nu [2l[ . 

Simulations have been performed at Pr = 1 and Ra = 
{10 5 , 10 6 , 10 7 , 10 8 } for an aspect-ratio T = L/H = 4 do- 
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FIG. 2: (Color online) Snapshot from one of the direct nu- 
merical simulations at Ra = 10 6 and Pr = 1. Shown is an 
iso-surface of temperature, colored by the kinetic energy. 

main. The grid resolution and other relevant information 
is given in Table HI The Reynolds number Re has been 
obtained from the peak of u'u' and Re T = u T H/is, with 

V ' T = \J v 'iz^ 1 ^ 2 a ^ ^ nc wau - Here, k represents the tur- 
bulent kinetic energy, which may not be the most ideal 
approximation of the shear velocity; normally the mean 
horizontal velocity is used. However, from the 'classical' 
(ensemble-average) point of view, there is no mean wind 
so that the only available data is from fluctuations. 

The grid clustering in the near-wall region has been 
chosen such that on average 8 cells were present in the 
thermal boundary layer. The kinetic boundary layer 
which is thicker than the thermal boundary layer at 
Pr = 1, contained about 16 cells on average. A snap- 
shot of one of the simulations at Ra = 10 6 clearly shows 
the unstable sheet-like plumes emerging from the bound- 
ary layers (Fig. [5]). Ten independent simulations with 
slightly perturbed initial conditions have been performed 
for all but the highest Ra, as the computational demands 
were too high. At Ra = 10 s on the 640 2 x 320 grid, one 
convective turnover time took 2500 hour on one SGI Ori- 
gin 3800 processor and even with 128 processors this is 
20 wall-clock hours per turn-over time. 

III. CLASSICAL RESULTS 

Instantaneous cross-sections of the temperature field 
are shown in Fig. [3] at Ra = 10 s . The dynamic be- 
havior can be viewed in the online animations [39j ] The 
vertical (x-z) cross-section of the temperature field (Fig. 
GJi) clearly shows the spatial segregation of hot areas 
where upward thermals dominate and cool areas where 
the downward thermals dominate. Fig. [3] shows a hori- 
zontal (x-y) cross-section of the temperature field at the 
edge of the thermal boundary layer. The boundary layer 
is a network of sheet-like plumes, which is coarse where 
the average flow is downward and dense where it is up- 
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(b) 

FIG. 3: (Color and movies online) Cross-sections of the tem- 
perature field at Ra = 10 s and Pr = 1. a) An x — z cross- 
section, b) An x — y cross-section at the edge of the bottom 
thermal boundary layer. The online movies are accelerated 5 
times, and blue and red represent low and high temperatures, 
respectively. 



wards. The sheets are formed by impingement of cold 
plumes onto the plate, as the hot fluid in the boundary 
layer is pushed away. These hot sheets move towards 
the region with ascending flow, where they seem to form 
an ever-contracting network of plumes. Where the net- 
work is dense, the plumes detach and the average flow is 
upward. 

Fig. shows the behavior of Nu as a function of 
Ra. This result is in good agreement with the relation 
Nu = 0.186 Ra 0276 , obtained byDNS with a similar do- 
main and boundary conditions [2l| , and with the classical 
wide-aspect ratio experiments of Chu and Goldstein [35| . 
The scaling of Re as a function of Ra (Fig. 3Jd), where 
Re is obtained from the maximum of u'u' , has a best- fit 
scaling as Re u = 0.17 Ra ' 49 . This is close to Re cx Ra 1//2 
which corresponds to a Reynolds number based on the 
free-fall velocity Uf = \/]3gKQTT. Note that the above 
scaling for Re is not presumed to describe asymptotic 
behavior, which cannot be expected in the range of Ra 
we consider. Instead it should be treated as a best-fit 
relation or local exponent. 
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TABLE I: Simulation details 



Ra 


grid 


At/t* x 10 3 


T/t* 


#sims 


Nu 


Re 


Re T 


1.15 x 10 s 


128 2 x 64 


1.13 


68 


10 


4.5 


54 


32 


1.0 x 10 6 


192 2 x 128 


0.57 


20 


10 


8.3 


157 


70 


1.0 x 10 7 


256 2 x 256 


0.45 


20 


10 


16.1 


458 


160 


1.0 x 10 s 


640 2 x 320 


0.11 


5 


1 


31.1 


1499 


210 



100 
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FIG. 4: a) Ra-Nu scaling for present DNS simulations, b) Re — Ra scaling for present DNS for Re based on horizontal squared 
mean fluctuations, along with a best-fit powerlaw. 



IV. WIND-DECOMPOSED RESULTS 
A. The wind structure 

In order to obtain the realizations for the symmetry- 
accounting ensemble-averaging, the complete three- 
dimensional fields for Ui, have been stored twice ev- 
ery convective turnover time, thereby ensuring that the 
fields are approximately independent. Furthermore, by 
performing different simulations at identical Ra with dif- 
ferent initial conditions, a real ensemble averaging was 
carried out. The realizations have been selected such 
that the wind structure has fully developed [l?], IH, H(| , 
so that the criterion for symmetry-accounting ensemble- 
averaging was satisfied. Over all ten simulations this 
resulted in approximately 400 independent realizations, 
which were then processed using symmetry-accounted 
ensemble-averaging, described in section lITBl 

The result of the averaging is shown in Fig. O for the 
simulations at Ra = 10 6 . Instead of a one-dimensional 
temperature profile &{z), a fully three-dimensional tem- 
perature field 0(x,y, z) is obtained of which an iso- 
surface is shown, clearly revealing the wind structure. 
These are the fingerprints of the role-like behavior of the 
wind structure. This is even better visible when mak- 



ing a slice through the hydrodynamic boundary layer 
(Fig. [5J3) . The contour lines arc of relative temperature 
r , which is the deviation from the plane-averaged tem- 
perature (Q)a(z), defined as Q r {x,y,z) = 0(x,y, z) — 
(&)a{z). The relative temperature 0,. is closely re- 
lated to the height-averaged temperature (&)h when 
(0)y = 0, as (6 r )ff = (0)h- The relative tempera- 
ture r is an indicator for where the fluid is rising and 
falling, as can be seen from the streamlines of the hor- 
izontal components u, v. Figure [5^ shows a side-view 
of the average field after averaging over the y-direction. 
Again, the iso-contours are of relative temperature r . 
Clearly visible in the figure is the projection of the two 
rolls onto the side view. Note that the periodic boundary 
conditions rule out the one-roll wind structures that are 
common for small-aspect ratio cells because of continuity 
arguments. 

In Fig. [6l the correlation of the height-averaged tem- 
perature (0) h with the wind structure (Q)h is shown 
as a function of time for the ten independent simula- 
tions at Ra = 10 6 . As will be recalled this is the match- 
ing criterion for the symmetry-accounted average, so the 
correlation with (O)h is an indication of how appropri- 
ate the method is, and also for the strength of the wind 
structure. It can be seen that on average, the correla- 
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FIG. 5: Results after symmetry-accounted ensemble-averaging at Ra = 10 6 and Pr = 1. a) (color online) 3D iso-surface 
of temperature, colored by the kinetic energy; b) Plane-cut in the hydrodynamic boundary layer, iso-contours of relative 
temperature Q r and streamlines of the horizontal velocity components; c) Result after averaging over the y-direction (top to 
bottom in Fig [5p). 



Ill 


1 1 


















III 


i i 



5 10 15 20 25 30 35 

t/t* 



FIG. 6: Correlation with the wind structure for the 10 simu- 
lations at Ra = 10 6 and Pr = 1 as a function of time. 

tion C with the wind structure is quite good, fluctuating 
between 0.5 — 0.85 for all simulations. 



B. Plane-averaged profiles of kinetic energy 

Plane-averaged profiles of kinetic energy k(z) = 
( u 'i u 'i)A and its components are shown in Fig. [7J Only 
one of the horizontal components is shown due to homo- 
geneity. The classical statistics (Fig.|7Ji,c,e) only differen- 



tiate between the horizontal and vertical fluctuations as 
the average velocity ul = 0. For this reason all variance 
of the wind structure is transferred to the fluctuations. 
From Figs. [7^,c,e one gets an image in which near the 
bottom wall variance of u'u' is created due to the ac- 
tion of the plumes impinging on and ejecting from the 
boundary layers. The interpretation from the symmetry- 
accounted profiles (Fig. [7b,d,f) is completely different. 
Here one sees that the maxima in u'u' = (uu)a + (u'u') a 
are primarily caused by the wind. The fluctuations, rep- 
resenting the action of the plumes, are nearly uniformly 
distributed in the bulk of the flow, and only a slight in- 
crease is visible near the boundary layers. The profiles 
of Fig. [7J scale nearly perfectly with the squared free-fall 
velocity Uj = /3gAQH for all three Ra numbers. Note 

that the plane-averaged momentum flux (w'u')a is not 
included in Figs. [7Ji-f, as this term is zero due to the 
symmetry of the wind structure. 

C. How does the wind affect the heat transport? 

Identifying the wind changes the decomposition of the 
vertical heat-flux. Using jJJ), ([8]) and the fact that (w)a = 
throughout the domain, we can rewrite (j4|) as: 

Nu = [(we) A - (w~t@>) A - K d z (e) A ) ■ (12) 
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FIG. 7: Plane-averaged profiles of kinetic energy. Shown are the classical profiles (a,c,e) and the symmetry-accounted profiles 
(b,d,f). a,b) Ra = 1.15 x 10 5 ; c,d) Ra = 1.0 x 10 6 ; e,f) Ra = 1.0 x 10 7 . 



As Nu is constant, only the distribution of the three terms 
on the right-hand side can change as a function of z. 
This is shown for the simulation at Ra = 10 6 in Fig. 
Diffusive transport dominates in the boundary layer, 

where the heat is transfered to the fluctuations w'& by 
entrainment/detrainment. In the bulk, about 30 percent 
of the heat is transported by the wind. Here we note 
that a simple model using sheet plume parameters (37j 
also yields that 30% of the heat is transported by the 
mean flow at Ra = 10 6 . 

Where the wind impinges on the wall, the boundary 



layer will be compressed and the local Nu will increase. 
Similarly, the local Nu will decrease in detachment zones. 
This effect is demonstrated in Fig. [5Jd where Nu as a 
function of x for the y-averaged wind structure (Fig. [5J;) 
is shown for the top- and bottom wall. Note that the 
spatial variations in the wall heat-flux are generated en- 
tirely by the wind structure since Nu(x, y) = — -£@d z Q 
at z = ±H/2. It can be imagined that spatial variations 
in Nu indicate significant horizontal heat-fluxes as well. 
Indeed, this is the case and this point will be addressed 
below. 
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FIG. 8: Balance ol heat-fluxes at Ra = 10 and Pr = 1. a) 
wind-decomposed heat-fluxes, b) Nu as a function of x and 
averaged over y at the top and bottom wall. 
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FIG. 9: Horizontal heat-fluxes are larger than the vertical 
and dominate deep in the thermal boundary layer, a) Vectors 
of the total convective heat-flux {uiO) y + (uj©') B and iso- 
contours of relative temperature at Ra = 10 6 . b) zooming in 
onto the boundary layer. 



The average horizontal heat-fluxes (uiQ)a and (u'jQ^a 
for i = {1,2} are zero by definition due to the absence 
of a forcing in the horizontal directions. However, as can 
be seen in Fig. [9^, where the total convective heat-flux 
(averaged over the y-direction) is shown in flux-vectors 

(u0 + u'Q',w& + w'Q'), the horizontal heat-fluxes are 
significant, especially very close to the walls. The heat 
transport is in the same direction at the top and bottom 
plates, and is directed to the relatively hot region where 
the flow is upward on average. 

Due to the anti-symmetry of uQ and u'& (Fig. (5Jt), 
their plane-average vanishes. Hence, (56)a and (u'Q')a 
cannot be used as an indicator for the strength of the 
horizontal heat-flux. However, the spatial standard devi- 
ation <7~g and ct^tq, are good indicators, with ax defined 
as 

a x = y/((X-(X) A )*) A . (13) 

The spatial standard deviations (Fig.O)) emphasize how 
close near the wall this heat is transported: The peak of 
the horizontal heat transfer lies deep inside the thermal 
boundary layer. This peak originates purely from the in- 
teraction of the mean wind and mean temperature field 
as u0. Horizontal heat-fluxes even exceed the average 



vertical heat-fluxes. These findings emphasize the im- 
portance of understanding the boundary layer structure 
and its dynamics. 

The error bars around the total heat-flux denotej^he 
spatial variations in the total vertical heat-flux wQ + 
w'Q' — nd z Q. An interesting aspect is that these varia- 
tions are large near the walls (due to the spatial variations 
in Nu, see Fig. [8)d), decrease and go to a minimum at 
z = Ae, after which the variance increases again due to 
the turbulent fluctuations. This suggests that the ther- 
mal boundary acts MS c\ redistributor of heat. 

The horizontal heat fluxes become larger as Ra in- 
creases, as shown by the characteristic heat fluxes nor- 
malized by Nu in Fig. [TTJ] Shown are the characteristic 
heat-flux due to the interaction of mean wind and tem- 
perature a~Q and turbulent heat flux ct^Tq, ■ Although 
the fluctuations ct^tq, grow in strength relative to Nu 
as Ra increases, their magnitude is still quite small at 
Ra = 10 7 . In contrast, the heat-flux due to the wind <t~q 
is nearly a factor 3 larger than the vertical heat flux at 
Ra = 10 7 ! The horizontal heat fluxes are central to the 
mechanism driving the wind, as is discussed below and 
in section fVl 
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FIG. 10: The peak of spatial standard deviation of the hori- 
zontal heat-fluxes normalized the heat-flux at the wall (<j))\ w 
as a function of Ra, showing the increasing wind-induced hor- 
izontal heat transport. 



TABLE II: Budget terms for momentum and heat equation. 
A V V B K 

d t ui — —djUjUi +vd 2 jUi —dip +/3gQ8i3 —dju'jU^ 



D. A wind feedback cycle 

In this section we study the momentum- and heat- 
balances term by term ( Table HI)). As the wind structure 
is statistically in a steady state, the balance is purely 
a function of space as A + V + V + B + TZ = 0. Sim- 
ilar to Fig [St, the budget terms have been averaged 
over the y-direction for convenience of presentation. Sev- 
eral checks were done to ensure that the y-averaged mo- 
mentum budgets are also representative for the three- 
dimensional field. 

In Fig. [TTIfour vertical sections arc shown, at the loca- 
tion of maximum upward motion (Fig. lllb ,). at 1/3 of the 
cycle (Fig. [TTb). at 2/3 of the cycle (Fig. Gib) and at the 
maximum downward motion (Fig. lllH ). Note that this is 
only half of the flow field; the other half does not provide 
new information due to symmetry. In the description it is 
sufficient to focus on the top wall only, as the top profiles 
from Fig. [TTa can be mapped onto the bottom profiles 
from Fig. Illt l by elementary symmetry operations, and 
the same holds for Fig. [TTb and Fig. [Tib . Focusing on 
the region where the flow is upward (Fig. \TTh). the forces 
of the horizontal momentum equation are nearly zero. In 
the vertical momentum equation, the buoyancy term B 
is balanced by the vertical pressure gradient V and the 
Reynolds stress 1Z. In this region, the average tempera- 
ture is positive, resulting in a positive buoyancy forcing B 
over nearly the entire vertical. The vertical pressure gra- 



dient is negative with a negative peak near the top plate 
which reflects the resulting pressure build-up due to the 
impinging plumes. The Reynolds stresses 7Z, dominated 
by the term —d z w'w', are slightly stronger on the top 
plate than on the bottom plate. This is an indication 
that on average, plume impingement is a more violent 
process than plume detachment. In the budget for tem- 
perature, the balance is primarily between diffusion T>, 
gradients in the turbulent heat-flux 1Z, with a small con- 
tribution due to the local acceleration of the mean flow 
field A. The forcing is stronger at the top plate, due to 
the impingement of the wind and the plumes. The local 
Nusselt number Nut is maximal at this position (see also 
Fig. [5t>). Note that Nu t is related to the integral of the 
thermal diffusive term T> on the top boundary layer. As 
the area under T> at the top-wall is larger than the area 
under V at the bottom wall, it follows that Nut > Nub, 
which is consistent with Fig. [5J 

Following the flow along the top plate, the horizontal 
momentum budget of Fig. [TTb shows a strong positive 
horizontal pressure gradient V , which is balanced by dif- 
fusion T> close to the wall, Reynolds stresses 1Z and iner- 
tial terms A a bit further away. The horizontal pressure 
gradient V is positive over the upper two thirds of the 
vertical. The interesting small peak in 1Z very near the 
wall will be discussed in more details in the accompany- 
ing paper pp| . which focuses on the boundary layers. In 
the vertical momentum equations, the situation is similar 
to that of Fig. [TTb . with the exception that the buoyancy 
force has become less positive. For the temperature bud- 
get, Nu t is lower at this point here (Fig. [5b), making 
thermal diffusion T> weaker. 

A bit further downstream (Fig. [TTb). the horizontal 
momentum budgets indicate that the pressure gradient 
is still positive but has decreased in strength. As the 
flow has started to decelerate, the inertial force A has a 
positive contribution. Close to the wall, diffusion T> is 
braking the fluid, and a bit further away the fluctuations 
1Z. As far as the temperature budget is concerned, Nu t 
has decreased even more. The budgets when the flow 
comes to a halt and starts its descent down are shown 
in Fig. IllH . In the vertical momentum equation, the 
buoyant forcing has become negative over nearly the en- 
tire vertical, which is balanced by the vertical pressure 
gradient V and the Reynolds stress term 7Z. As Nu* is 
at a minimum at this position, thermal diffusion is rela- 
tively small here, and the advective part A has become 
negligible. 

Concluding, the mean momentum and temperature 
budgets show that the wind is driven by pressure gra- 
dients. These pressure gradients are generated as the 
result of spatial buoyancy differences caused by spatial 
temperature differences. This finding is in line with the 
study by Burr et al. [25[ , despite the absence of sidewalls. 
The pressure gradient can be estimated by integrating 
the vertical momentum equation, as will be shown in the 
next section. 

Using Fig.[JTJwc can identify a detailed feedback mcch- 
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FIG. 11: Momentum and temperature budgets as a function of z for Ra = 10 6 and Pr = 1. a) upward motion; b) 1/3 of the 
way; c) 2/3 of the way; d) downward motion. Note that only half of the wind-structure shown in the center picture (see Fig. 
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FIG. 13: Sketch of the wind structure and 9 locations A-I. 



anism sustaining the wind. The buoyancy force creates 
a pressure increase (decrease) on the top wall where the 
flow is positively (negatively) buoyant. This generates 
horizontal pressure gradients at the top- and bottom 
walls that drive a mean flow which transports a relatively 
large amount of heat through the bottom layers (section 
IIV C[) . The net transport of heat towards the region with 
ascending flow causes spatial temperature gradients (Fig. 
[9]). Finally, these spatial temperature differences induce 
spatial gradients in the buoyancy which completes the 
feedback cycle. A schematic diagram of this process is 
shown in Fig. [T2] 



V. A SIMPLE MODEL FOR THE WIND 

A. A short derivation 

Based on the feedback mechanism deduced in the pre- 
vious section, a simple mathematical model can be con- 
structed, by averaging the two-dimensional momentum 
and temperature equations over appropriate regions of 
space. A sketch of a typical wind structure is shown in 
Figure [131 with 9 locations A-I which will be used to 
identify specific areas. A generic averaging operator (■), 
which averages both over lines and areas, is defined as 



(X)ci=^Jxdz, 
(X)acig 



Xdxdz, 



and so on. Here, L w represents the size of a roll (Fig. 
[T5)) . As there is a slight clash of variable names (with the 
height H) , it should be understood that the locations A-I 
will only be used as subscripts in the averaging operator. 

The model has two main variables, the mean wind ve- 
locity U w and the mean temperature amplitude Q w . The 
mean wind velocity U w is defined as 



{UjACFD 



HL, 



udxdz. 



J w JJACFD 

The mean temperature <d w is defined as 



(6) 



BCIH 



HL, 



Qdxdz, 



(14) 



(15) 



BCIH 



which represents the wind-induced temperature ampli- 
tude. 

Averaging the two-dimensional horizontal momentum 
equation over the area ACFD and the temperature- 
equation over the area BCID results in 



dU w _ __^(w'u') DF 



CF 



dt 



H 



)ad _ „ {d z u)AC 



H 



(16) 

de w 2(uO) BH , 2(^e') BH (d z e) H i-(d z e) B c 

~ ; 1 h K 



dt 



H 



(17) 



A technical discussion about the steps leading to 
([17]) can be found in appendix lAl 

In the horizontal momentum equation (|16p , we see that 
U w is driven by a yet unspecified pressure gradient, and is 
subject to a wall shear stress and a turbulent shear stress 
in the bulk (see Fig. 14(a)). Both the wall shear stress 
and the turbulent stress tend to decelerate the wind. In 
the heat equation (|17|). the temperature amplitude Q w 
is driven by the large horizontal heat flux (uQ) bh in the 
boundary layer, which was identified in section lTV CI The 
term (u'Q')bh is a horizontal turbulent heat-flux, which 
tends to decrease temperature differences by turbulent 
mixing. The last term in (|17|) represents the heat-flux 
through the bottom and top wall. If Q w is positive, the 
heat-flux on the top-wall will be larger than the heat-flux 
on the bottom wall. Hence, this term effectively removes 
heat from the co ntrol v olume. A sketch of the heat-fluxes 
is shown in Fig. |14(b)| 

The average pressure gradient can be estimated with 
the help of the vertical momentum equation. Averaging 
the vertical momentum equation over CI, which is the 
streamline connecting the bottom to the top wall, results 



d(w) 



ci 



dt 



Pg(Q) 



ci 



c 



H 



(18) 



w J JACIG 



Thus, the average vertical acceleration over CI depends 
on the average temperature and the pressure difference 
between the top and the bottom wall. Because of the 
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FIG. 14: a) Dominant forces on the wind structure. A pres- 
sure gradient drives the wind, while the wall-shear stress and 
turbulent shear stress in the bulk provide friction; b) Heat 
fluxes due to the wind structure. The heat-flux uQ creates 
spatial temperature differences, while the heat-flux at the top 
and bottom wall and the turbulent heat flux in the bulk de- 
stroy temperature differences. 



point symmetry around E (Fig. I13p . the pressure (p)i 
is equal to (p)a, which means that (|18p provides infor- 
mation about the mean pressure gradient on the bottom 
wall. Invoking continuity and approximating the pres- 
sure gradient as a linear function of z (see appendix |A")) 
yields 



ICF 



)ad H 2 dU u 



2L 2 „ dt 



fan 

2L,„ 



(19) 



This is one of the central results of this paper, as ([19] 
provides an explicit coupling between U w and Q w . 

Substituting (fl"9j) into (fl6|) yields the unclosed equa 
tions governing the wind structure: 



dUw 
dt 



2L 2 



2L 2 + H 2 \2L 



H 




d@ w 2{uQ) BH , 2{u'Q')bh , (d z G) H i- (d z G) BC 

- 1 h K 



dt 



H 



(21) 




FIG. 15: Generation of the horizontal heat-flux u0, which 
generates spatial temperature differences. 



B. Parameterization, turbulence closure and 
dimensionless formulation 

The viscous momentum and diffusive heat fluxes at the 
walls in ([21?]) . (|2"Tj) can be related to U w , <d w and Ae by 

v{d z u) A c ~ \ U ™\ Uw 



K(d z Q) H i ~ k- 

n(d z Q)BC ~ K 



-Ae/2-9^ 

e w - Ae/2 

A^ ' 



The wall shear stress v(d z u)AC is expressed simply in 
terms of the friction factor Cj [38j . The temperature 
gradient at the top wall (d z Q)Hi can be estimated by 
(— A0/2 — Q w )/Xq, as variations in Ae are negligible to 
first order. The temperature gradient at the bottom wall 
is approximated similarly. The mean horizontal heat-flux 
(uQ) bh which drives the flow (section IIVC|) , is approx- 
imated by 



(50) 



BH 



H 



The horizontal heat-flux occurs mainly in the thermal 
boundary layers (Fig. I15p , where the temperature is ap- 
proximately A0/2 and the typical velocity is U w . Hence, 
uQ f=a U w AQ/2, and accounting for the two bound- 
ary layer contributions, the average horizontal heatflux 
(uQ) bh is approximated as above. 

The only terms which require closure at this point are 
the turbulent momentum and heat flux, (w'u')df an d 
(u'Q')bh respectively. The bulk is well-mixed, as can 
be judged from the nearly constant temperature and the 
linearly varying velocity as a function of z in the bulk. 
Therefore, a simple closure with the gradient-diffusion 
hypothesis is appropriate for the turbulent fluxes 



(u'& 



(w'u')df = -v T d z u : 
bh ~ -Kr{d x &)BH = 



2f4, 
H 

Pr T L w 



(22) 



(23) 
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where vt and Pr-r = vt/kt are the eddy viscosity and 
turbulent Prandtl number, respectively. To relate i>t to 
mean flow properties, we use the Prandtl mixing length 
hypothesis, which results in 



v T = a£ 2 \d z u\ w aH : 



H 



a\U w \H 



(24) 



Here a is a free parameter which controls the mixing. 

Using the approximations above, the equations for the 
wind structure are given by 



dUu 
dt 



2L 



2LI + H 2 \ 2L W 



4a + C f 
H 



UJU U 



dQ w 2A e A6 4a\Uy\H 2 



dt 



(25) 
(26) 



Introducing dimensionless variables U w = U w /Uf, 
® w = Q w /A@, i = tUf/H, where U f is the free-fall 
velocity Uf = ^(3gHAQ, results in 



dU u . 
di 

dB^ 
di 



2L 



/ 1 



2i2+i V2L, 



■e, 



(4a + C f ) 



(27) 



2A e 



4o 



e, 



A e Re/Pr 



(28) 



Here, L ro = L w /H and Ae = X&/H are the normalised 
roll size, kinetic and thermal boundary layer thickness. 
Re/ is the Reynolds number based upon Uf. 

The wind model (|27( I28| comprises two nonlinear cou- 
pled ordinary differential equations in U w and Q w . The 
model contains seven parameters, L W1 Cf, Ae, ct, Re/, Pr 
and Pry. However, Ae = Ae(Ra, Pr), Cf = C/(Ra, Pr) 
and Re/ = Ra 1 ^ 2 Pr -1 ^ 2 . Therefore, the model can be 
expressed the parameters Ra, Pr, L w , a and Pry com- 
plemented by the functions for X u and Cf. Only Pry and 
a can be used to calibrate the model, which will be done 
based on the simulations at Ra = 10 6 in the next section. 



C. Results 

In this section the model will be compared to the DNS 
results. As a baseline test, the wind model (|27l |2"5)) 
should be able to predict the trends in wind speed U w 
and temperature amplitude Q w as a function of Ra. In 
this study, we close L w , C / and Ae empirically with our 
DNS results. In particular, we use Pr = 1, L w = 2\/2, 
Cf = A T Rai- and A e = A e Ra 7e . The best-fit coeffi- 
cients for Cf and Ae based on the current simulations 
are A T = 36, A e = 2.33, j T = -0.30 and 7e = -0.27. 

The turbulence parameters a and Pry will be cali- 
brated using the turbulent fluxes and wind and temper- 
ature amplitude for the simulation at Ra = 10 6 . By 
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FIG. 16: Phase-space of wind model at Ra = 10 7 . The fixed 
points are denoted by circles and the black line is the separa- 
trix. 



calculating v t and kt with (|22j) and (|23|) it follows that 
Pry ~ 0.85, in reasonable agreement with the generally 
accepted Pit « 0.9 for shear flows (38|. The mixing pa- 
rameter a can be calculated from (|24j) . which results in 
a « 0.6. It is noted that a and Prr are not parameters 
in the strict sense, as the DNS results indicate they have 
a weak dependence on Ra. 

The phase-space of (|27l |28|) at Ra = 10 7 is shown in 
Fig. [T7JJ There are three fixed points in the domain, of 
which the one at (0, 0) is a saddle node. The two other 
fixed points are attractors. Thus, if there is no wind ini- 
tially, any small perturbation caused by turbulent fluctu- 
ations will cause the system to settle in a wind structure 
with either U w > or U w < 0. The tendency of Rayleigh- 
Benard systems to establish a wind structure can thus be 
explained by the positive feedback created by wind ad- 
verting large amounts of heat and the resulting buoyancy 
differences which drive a mean flow. The amplitude of the 
wind is the result of the interaction between the destabil- 
ising mechanism mentioned above and the mixing due to 
turbulence which reduces gradients. Note that the model 
cannot describe wind reversals @, [H, [H| , by the absence 
of dynamic fluctuations; both nonzero fixed points are 
stable. The limitations of the model will be discussed in 
more details in the concluding remarks (section IVlj) . 



As the system is invariant under U w — > —U w , Q w — > 
Q w it suffices to study the positive fixed point of (|2"T1 
which is located at 




e -— u 2 

ai 



(29) 
(30) 
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FIG. 17: Behavior of the model (eqns ([29 
compared to the DNS data (diamonds), a) 
of Ra; b) O m as a function of Ra. 



where 



, (J30J) , solid line) 
U w as a function 



2L, 



2A f 



-, a 2 = 4a + Cf, 



4a 



A e Re/Pr 



Shown in Fig. 17(a) and 17(b) are the trends of U w and 



Q w as a function of Ra, compared with the DNS results 
(diamonds). The model slightly underpredicts U w , but 
the temperature amplitude Q w is predicted well. More 
importantly, the model seems to capture the decreasing 
trend of <d w properly, as well as very weak Ra dependence 
of U w . Given its simplicity, the model is in fair agreement 
with the simulations. 

From (f2T))) . it follows that as Ra increases, Cf becomes 
negligible relative to the mixing parameter a. For the 
simulation at Ra = 10 7 , Cf w 0.17 while 4a = 2.4. 
Hence, the friction term Cf + 4a is dominated by the 
turbulence in the bulk. As Cf is a decreasing function 
of Ra, this effect becomes stronger as Ra increases. This 
indicates that wall friction has a negligible influence on 
the wind velocity for Ra sufficiently high. 

The asymptotic scaling of U w for Ra — » oo can be 



established by studying the scaling of the coefficients of 



b 2 
b\b 2 
a 2 b\ 



A e L 2 w Pr T 
2aA e 
Act 



Ra 



-(l/2+ 7e ) 



Pr 



-1/2 



4« + Cf Li Pr 



A% Ra l+3 7 ep r 



Assuming that the scaling exponent for Ae remains above 
7e = —1/3, Ra 1+37e — > oo, by which asymptotic scaling 



7e 
ofU 



is 



U w oc Ra 7e/2 



(31) 



The wind Reynolds number Re^ = UfHv 1 U w /Uf = 

U w Ra 1/2 Pr- 1/2 , so that Rc w oc Ra (1+7e)/2 . Based on the 

exponent from the simulations (7© = —0.27) it follows 

that Re„, oc Ra ' 37 in the asymptotic limit. As Ke F oc 
,0.44 



-7e)/3 



oc Ra ' 44 (where we used that Ae oc Nu 1 ), 



Ra^ 1 - 

Re^, oc Ra 3 ' suggests that the wind becomes progres- 
sively weaker relative to the fluctuations as Ra increases. 
Naturally one should not assign too much value to the 
exact exponent, as it critically depends on the Ra de- 
pendence of Ae. Nevertheless, the flux term generating 
temperature differences uQ depends critically on Ae- If 
Ae is a decreasing function of Ra, so will Q w and U w . 



VI. CONCLUDING REMARKS 

The aim of this study has been to clarify the processes 
responsible for the wind amplitude. Direct numerical 
simulation was performed at Ra = {10 5 , 10 6 , 10 7 , 10 8 } 
and Pr = 1 for an T = 4 aspect ratio domain with peri- 
odic lateral boundary conditions. For all but the highest 
Ra, 10 independent simulations were carried out, result- 
ing in approximately 400 independent realizations per 
Ra. The wind structure was extracted by accounting 
for symmetries, i.e. using the translational invariance of 
the system to align realizations before averaging them. 
In this way, wind could be distinguished from fluctua- 
tions for a domain with periodic sidewalls. It was found 
that the characteristic peak in the kinetic-energy pro- 
file by which the boundary layer thickness is defined, is 
nearly entirely due to the wind and the turbulent fluctu- 
ations (u'u')a arc distributed uniformly outside the ther- 
mal boundary layer. Deep inside the thermal boundary 
layers, the wind structure is responsible for large hori- 
zontal heat-fluxes, transporting heat towards the region 
of upward flow, through the terms uQ and vQ. These 
horizontal heat-fluxes are up to three times larger than 
the average Nusselt number at Ra = 10 7 , although the 
total amount of heat transported through the boundary 
layer decreases with Ra. This wind-generated horizontal 
heat-flux is central for the formation of a wind structure 
as it generates spatial temperature differences. As a re- 
sult of the temperature differences, pressure gradients arc 
generated which drive the wind. 
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A simple model of two coupled nonlinear ordinary dif- 
ferential equations was derived, which captures the es- 
sential processes governing the wind structure. The pri- 
mary variables are the wind velocity U w and the tem- 
perature amplitude Q w , while the Rayleigh number Ra, 
the Prandtl number Pr, wind roll size L w , friction fac- 
tor C/(Ra, Pr) and thermal boundary layer thickness 
Ae (Ra, Pr) are physical parameters. The turbulence in 
the bulk is described by a mixing coefficient a and a tur- 
bulent Prandtl number Pr^. DNS results were used to 
calibrate a and Pr^, and served as inspiration for the 
parameterisation. The model reproduces the Ra depen- 
dence of U w and Q w from the DNS, and the following 
conclusions follow from the wind model: 



APPENDIX A: DERIVATION OF WIND MODEL 

In this appendix we average the two-dimensional mo- 
mentum and temperature equations over specific control 
volumes in order to develop a theoretical model for the 
wind. The model has two variables, the wind velocity U w 
and the temperature amplitude Q w , which are defined in 
section [V] To identify different regions, various locations 
are denoted by A-I in Fig. [131 The wind roll size is de- 
noted by L w , and (•) is the generic averaging operator 
defined in section IVl 

1. Horizontal momentum equation 



• A wind structure is inevitable, as the fixed point 
corresponding to the absence of wind is an unstable 
saddle. The positive feedback responsible for this 
behavior is the interaction between the mean wind 
and the mean temperature, as described above. 

• The influence of the wall friction on the wind ve- 
locity is rather limited. At Ra = 10 7 , we find that 
Cf = 0.17, while 4a = 2.4, so that the turbulence 
in the bulk dominates the total friction Cf + 4a in 

Although the model gives interesting insights, it has 
a number of limitations. In the derivation it has been 
assumed that the domain was unbounded in the lateral 
directions, i.e. no sidewalls. As a result, the effect of fric- 
tion on the sidewalls has been omitted, which - once in- 
cluded - will enhance the friction experienced by the wind 
structure. Furthermore, the model was derived from the 
two-dimensional Reynolds-averaged Navier-Stokes equa- 
tions, which accounts only for the mean effects of the 
turbulence, thereby excluding long-term dynamical be- 
havior such as reversals and reorientations. However, no 
fundamental difficulties are expected to incorporate the 
missing physics described above. 

In the accompanying paper (26[, we focus on the 
boundary layers. Using the wind model developed in 
this paper, we derive new scaling laws for X u and Cf. 
For the wind model, this implies that Ae is the only free 
parameter in the wind model. Furthermore, we discuss 
in detail the issue whether or not the boundary layers 
should be regarded laminar or turbulent. 



Acknowledgments 

This work is part of the research programme of the 
Stichting voor Fundamcntccl Onderzoek der Materie 
(FOM), which is financially supported by the Ned- 
crlandsc Organisatic voor Wctcnschappelijk Onderzoek 
(NWO) The computations were sponsored by the Sticht- 
ing Nationale Computerfacilitcitcn (NCF). 



The two-dimensional horizontal momentum equation 
is given by 

dtu = — d x uu — d z wu — d x u'u' — d z w'u' 

(Al) 

-d x p + v{dlU + dlu). 

This equation will be averaged over the area ACFD, 
which results in 

=o =o =o =o 

dU w (uu) C F ~ (uu)ad {wu)df - {wu)ac 



dt 



H/2 



= 



=o 



(u'u')cf ~ (u'u')ad (w'u') df - (w'u')ac 



)CF - (P)AD 

RS0 



H/2 



(d x u) C F - (d x u) A D (d z u) DF -(d z u)Ac 
L w V H/2 

Due to the choice of the control volume, many terms are 
zero (indicated by — above them) . Other terms can be 
neglected (indicated by 0). The three viscous terms 
are neglected as they are very small compared to the 
wall friction term. The average horizontal fluctuations 
on the interface CF and AD will be approximately of the 
same strength, so that these terms cancel out. Hence, 
the horizontal momentum equation simplifies to 



dU u: 
~dt 



= -2 



[w'u') DF 
H 



/CF 



)ad _ 2 (d z U)AC 



H 



(A2) 



2. Temperature equation 

The temperature equation is given by 

d t Q = - d x uO - d z wQ - a^u 7 ©' - d z u/G' 
+ K{d 2 x ® + d 2 z Q). 



(A3) 
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This equation is averaged over the area BCIH (Fig. [T3|) . 
yielding 

=o =o =o 

dQ w (uQ)ci -(uQ)bh (wQ)hi - (wQ)bc 



dt 



L w /2 



H 



=0 =0 =0 

I^CI-{uW)bH (^Q')HI-(^e'}BC 



L w /2 



H 



(d x e)a-(d x e) B H , (d z &) H i - (d z Q) B c 

K ; \-K 



L w /2 



H 



Again, the choice of the control volume causes many 
terms to be zero (indicated by = 0), while other terms 
can be neglected (indicated by as 0). Here, the horizontal 
diffusive heat-fluxes can be neglected, because they are 
very small compared to the vertical diffusive heat-fluxes. 
The temperature equation is reduced to 

dQ w _ 2(uQ)bh , 2(u 7 6') bh , (d z Q) H i - {d z Q) B c 

7 1 7 r K 



dt 



H 



(A4) 



3. Continuity 

The continuity equation 

d x u + d z w = 0, (A5) 
is averaged over BCFE (Fig. [T3|) . which results in 



line CI (Fig. [13)) . As spatial derivatives in the unbounded 
directions are zero (see Fig. Ilia ), the vertical momentum 
equation reduces to 



d t w = [igQ — d z ww + d z w'w' + d z p + vd 2 w. (A7) 
Averaging over CI gives 



d(w) 



ci 



dt 



--/3g(0) 



ci 



)C 



H 



=0 =0 

(ww)i—(ww)c (w'w')i — (w'w')c 
H H 
=o =o 

{d z w)i - (d z w)c 
H 

It can be verified that d z w = at the bottom and top 
plate by substituting the no-slip boundary condition u = 
in the continuity equation. Hence, the average vertical 
momentum equation reduces to 



d(w) 



ci 



dt 



Pg(Q) 



ci 



)c 



H 



Due to symmetry, the pressure at A and I is identical. 
Hence, substituting (p)j = {p) a, estimating (w)ci ~ 
W w , (<d)ci ~ &w and using (|A6|) gives that the typi- 
cal pressure gradient at the bottom plate is given by 



(u)cf~(u)be , (w)ef-(w}bc 



L w /2 



H/2 



0. 



Estimating (u)be ~ U w and (w)ef ~ W w , with W w the 
mean vertical velocity, the continuity equation becomes 



Uu 



Ww 

H 



(A6) 



4. Vertical momentum equation 



(P)c - (p)a = H^dU^ _ PgH 
LI dt 



In Fig. [TTJ we can see that the pressure gradient is ap- 
proximately a linear function of z, by which the average 
pressure gradient can be estimated as 



(P)CF - (p)AD H 2 dU w /3gH 



L,, 



2LI dt 2L„ 



e, 



(A8) 



The unknown pressure gradient can be obtained by av- 
eraging the vertical momentum equation over the stream- 



Equations (|A2[) . (|A4[) and (|A8|) constitute the unclosed 
wind model. 
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